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The percolation threshold for flow or conduction through voids surrounding randomly placed 
spheres is rigorously calculated. With large scale Monte Carlo simulations, we give a rigorous 
continuum treatment to the geometry of the impenetrable spheres and the spaces between them. 
To properly exploit finite size scaling, we examine multiple systems of differing sizes, with suitable 
averaging over disorder, and extrapolate to the thermodynamic limit. An order parameter based 
on the statistical sampling of stochastically driven dynamical excursions and amenable to finite size 
scaling analysis is defined, calculated for various system sizes, and used to determine the critical 
volume fraction <f> c = 0.0317 ± 0.0004 and the correlation length exponent v — 0.92 ± 0.05. 

PACS numbers: 64.60.ah,61.43.-j,64.60.F-,61.43.Gt 
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I. INTRODUCTION 

Notwithstanding simple underlying Physics, percola- 
tion transition are bona fide second order phase transi- 
tions, with all of usual the hallmarks of singular behavior 
at a critical poinfci. Percolation transitions, phase tran- 
sitions mediated by a disordering influence such as the 
random removal of sites or bonds on a regular lattice ge- 
ometry, are amenable to numerical study in the context 
of Monte Carlo calculations used to treat disorder. How- 
ever, while the critical behavior in discrete systems may 
be characterized to a high level of precision in this man- 
ner, there are salient examples of percolation phenomena 
which cannot be reduced to discrete lattices. 

In realistic porous media, the flow of liquid on a macro- 
scopic basis typically entails the passage of liquid around 
irregularly shaped granules which comprise the material 
and represent barriers to fluid flow. For the movement of 
fluid or charge in porous materials, there is a dichotomy 
for volume concentrations of barrier particles above and 
below a critical concentration p c . If inclusions are suffi- 
ciently sparse, voids join together in a macroscopic con- 
nected network spanning the entire system, allowing for 
the transmission of fluid at the bulk level. On the other 
hand, for large enough concentrations of randomly placed 
barrier particles, contiguous navigable voids are limited 
in size to a finite length scale £, barring the transmission 
of fluid at a system wide level; the critical density p c of 
inclusions which separates the two scenarios is thus of 
practical significance. 

For the sake of convenience, we rescale coordinates in 
such a way that the spherical inclusions are each of unit 
radius. Although in principle one may record the critical 
density p c where the percolation transition occurs, a more 
frequent practice in discussing percolation through voids 
is to report instead the excluded volume^, given by <p c — 

Determining if voids may be traversed at the macro- 
scopic level is a continuum percolation problem for which 
a topologically equivalent discrete counterpart is in gen- 



eral not readily accessible, a condition very distinct from 
that of a complementary percolation problem (with the 
permeability pattern reversed) where the flux occurs 
through networks formed of overlapping included par- 
ticles with interstitial regions impervious to fluid flow. 
The geometry of the randomly placed barriers, known 
a priori, facilitates the calculation of characteristics re- 
lated to percolation when the flow is through included 
particles instead of the interstitial spaces. On the other 
hand, the fact that the shapes of voids between barrier 
particles are only evident a posteori hampers a rigorous 
identification of networks formed by connected voids and 
enhances the challenge of finding the percolation thresh- 
old. Nevertheless, with the aid of a finite size scaling 
analysis, we determine not only the critical volume frac- 
tion (p c , but also the critical exponent v associated with 
the correlation length £. 



Previous theoretical studies include calculations which 
use an approximate discrete scheme and extrapolate to 
the continuum case&i. In a related effort, the discretiza- 
tion scheme has been applied to randomly placed ellip- 
soids 5 and oblate spheroids 6 . In particular cases, such as 
systems comprised of randomly placed spheres, discrete 
networks equivalent to void spaces^ may be constructed. 
Voronoi tessellations, and the appropriate generalizations 
have been applied^— to systems comprised of randomly 
located spheres. In a recent set of calculations, there is no 
discretization scheme and the simulation entails following 
tracer particles in an effectively infinite sized syste m 11 ! 12 ; 
the percolation threshold is identified by determining the 
density of spheres for which the tracers cease to move 
diffusively. In this work, we use a dynamical approach 
where systems of finite size are examined. In this way, we 
are able to apply standard machinery of finite size scaling 
analysis while providing rigorous treatment to the geo- 
metric intricacies of voids between barriers which serve 
as conduits for fluid or charge in the percolative regime. 
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II. CALCULATION OF THE EXCURSION 
ORDER PARAMETER 

An order parameter, an intensive thermodynamic vari- 
able finite when the phase is intact, and decreasing con- 
tinuously to zero at a phase boundary (e.g. magnetiza- 
tion per unit volume in the case of ferromagnetism dis- 
rupted by thermal fluctuations above the Curie Temper- 
ature and nonzero below it), may be pressed into ser- 
vice as a tool to locate second order phase transitions 
and thereby determine the phase diagram. In percola- 
tion phenomena where connected navigable clusters are 
readily identified, the "strength" of the spanning clus- 
ter, or the number of sites in the lattice it encompasses 
relative to the total number of sites in the system, is 
conveniently sampled in Monte Carlo simulations using 
techniques such as the Hoshcn-Kopelman algorithm for 
the identification and characterization of connected clus- 
ters^. 

In the scenario of interest, where percolation occurs 
through irregularly shaped interstices among randomly 
placed barrier particles of a prescribed shape, the connec- 
tivity of void volumes to adjacent navigable regions may 
be difficult to establish in more general situations (e.g. 
non-spherical inclusions such as tetrahedra, cubes, or oc- 
tahedra). To address this challenge, we use dynamical 
simulations to probe and determine the extent of voids. 
Although stochastically driven exploration is validated in 
the case of randomly placed spheres, the approach may 
be used in scenarios where variants of Voronoi tessella- 
tion are more difficult to apply. 

In the spirit of a Gedanken experiment, we envision 
applying a reflective coating to exposed surfaces of over- 
lapping spheres such that light striking the included par- 
ticles is reflected with no penetration of the incident rays. 
Whereas the transmission of light through the matrix of 
inclusions would signal that the system would admit the 
flow of fluid or electrical charge, a configuration which 
blocks the transmission of incident light also lacks a per- 
colating path. 

Dynamical trajectories are initiated with an unbiased 
choice for the point of origin in the interior of a void and 
exterior to all barrier particles. To select a void point, 
randomly chosen candidate locations are accepted only 
if the point is exterior to all neighboring spheres. Some 
economy is gained in partitioning the 3D system into 
a cubic grid, where the unit cell size is slightly greater 
than the radius r of the spheres. The subdivisions reduce 
the task of checking for penetration into nearby spheres 
to simply checking the cubic unit cell occupied by the 
prospective (randomly chosen) starting point as well as 
each of the adjacent cubes. To set the virtual ray in 
motion, a "velocity" vector v is randomly chosen in a 
spherically symmetric way to eliminate any directional 
bias, with v of unit length since "speed" has no physical 
or computational significance in the context of the order 
parameter calculation. 

The ray trajectory is xi — xq + vt with the start- 



ing point xo and the orientation of the velocity vec- 
tor v stochastically determined. The path of the ray 
must be traced until contact is made with the near- 
est sphere; with a barrier particle centered at x c , the 
condition for intersection with a sphere of radius r is 
\xi — x c \ = r. Standard geometric arguments lead 

to dcoii = —XdCosO ± yl — x 2 d sin 2 for the distance 
to a point of intersection, where Xd = \xq — x c \ and 
cos 6 = —v-Sd- Choosing the negative sign corresponds 
to the closer point of entry and the positive sign to the 
more distant point of departure as the ray propagates 
through and exits from the sphere. Finding the closest 
included particle intersected by the ray is tantamount to 
minimizing cJ co ii, the shortest distance to a sphere along 
the trajectory. 

To reduce computational overhead, the system is sub- 
divided into a cubic grid in the same fashion as used to 
determine if a randomly chosen point is located in a void; 
rays are propagated from one small unit cell to the next 
until a collision with a barrier particle is recorded. At the 
site of a collision x co n, to avoid penetration of the imper- 
meable particles comprising the matrix, the direction of 
motion must change in such a way that the incident ray 
is projected away from the interior of the barrier parti- 
cle. Specular reflection is implemented by reversing the 
component of the velocity in the direction of x C oii — x c , 
the radius vector extending from the sphere center to the 
site of the collision. 

Although mirror reflection encapsulates the physics of 
geometric optics appropriate to a ray of light moving in 
a matrix of silvered spheres, propagation by simple spec- 
ular reflection may be improved on, particularly in the 
case of neck-like voids where many collisions may be re- 
quired to traverse the narrow corridor if the ray is di- 
rected normal to the primary axis of the void. To im- 
prove the efficiency of navigation by roughly an order of 
magnitude, we incorporate moves where, instead of mir- 
ror reflection, the ray trajectory shifts by a right angle. 
In 3D, even with a perpendicular shift imposed, there re- 
mains a degree of freedom, which is selected at random to 
provide stochastic input at each collision. The balance of 
the moves is such that stochastically infused right angle 
moves outnumber specular moves two to one. 

To interrogate the system and determine whether a 
void is part of a spanning cluster, we calculate the ex- 
cursion parameter 6 = j^{Sf ma + 5}! ms + 5* ms ), structured 
to allow for the incorporation of statistical information 
gleaned from traversals along each of the three Cartesian 
axes. In calculating the rms deviations, special consider- 
ation must be given the periodicity of the system. When 
the density of randomly placed spheres is high enough to 
block percolation, where <j> < (j> c , the scale £ of contigu- 
ous void networks is in most cases small relative to the 
system size, and accordingly, i5 < 1. 

On the other hand, if inclusions are sparse enough 
that <f> > <fi c , system wide traversal becomes feasible. 
Moreover, with periodic boundary conditions imposed 
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the RMS deviations along one or more axes could in prin- 
ciple grow without bound as the ray propagates from one 
replica of the system to the next. In such situations, the 
result would be an unrealistically large order parameter, 
a circumstance which may be remedied by truncating 
5f mB and counterparts along other axes at the system 
size L. Thus, in the extreme case of a system devoid of 
inclusions, 8 would attain its maximum value of 1. The 
excursion quantity <5, which saturates at unity in the ab- 
sence of obstacles and falls to zero at the percolation 
threshold for bulk systems, serves as an order parameter 
giving a measure of the size of a randomly selected void. 
In the thermodynamic limit, the configuration averaged 
8 also gives the probability that a randomly selected void 
is a system spanning cavity. 

In a diffusive excursion, though the disorder and time 
averaged vector displacement vanishes, the root mean 

1/2 

square distance traveled scales as N ' a , where N co n is 
the number of ray interactions (and concomitant trajec- 
tory deflections) with barrier particles. Unless one is pre- 
cisely at the percolation threshold with RMS deviations 
from the trajectory starting point scaling as A^ oll , where 
e < 1/2 in the non-diffusive critical regim o 11 ! 12 , configu- 
ration averaged ray excursions will scale diffusively with 
the number of sphere interactions, if the system is suffi- 
ciently large. 

If the void is small in size, a ray which propagates in 
the cavity will probe its extent relatively quickly, yielding 
a finite 8 value. On the other hand, with a system span- 
ning void scaling as the system size L, the ray traversal 
time varies as L 2 on average. Hence, to avoid a spuri- 
ously low result for the percolation threshold (f> c , trajec- 
tories are propagated with collisions accumulating until 
N co i\ 3> vL 2 where r\ is an integer prefactor. By succes- 
sively doubling 77, a saturation point is eventually found 
(i.e. for 77 = 80, 000 and 77 = 160, 000) where statistically 
significant changes in critical quantities cease to occur. 

The calculation of the 8 parameter entails averaging 
over many realizations of disorder (i.e. at least 10 ) to 
suppress statistical fluctuations and permit the obser- 
vation of finite size scaling trends, which are exploited 
to calculate critical quantities such as <p c and the cor- 
relation length exponent v. To generate disorder con- 
figurations, we sample the grand canonical ensemble as 
described in the Appendix to determine the number of 
spherical barriers encompassed in the simulation volume. 
The Mersenne Twister algorithm, which generates high 
quality random numbers with a modest computational 
cost, is used for stochastic input. 



III. FINITE SIZE SCALING ANALYSIS 

Percolation transitions, whether in the context of con- 
tinuum or in discrete cases, are second order transitions 
with singular behavior in thermodynamic variables at 
the phase transition. Single parameter finite size scal- 
ing analysis posits a form L a f([(f> — </> c ]L 1 / !/ ) for ther- 
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FIG. 1: (Color online) The excursion order parameter 8 is 
plotted with respect to excluded volume 4> for various system 
sizes; the inset is a magnified view of the intersection of the 
8 curves. 



modynamic variables, a singular form anticipated for the 
excursion order parameter 5. Since 8 is normalized and 
thus dimensionless, a = 0, and the singular dependence 
reduces to 8 ~ f{[4>—(f> c \L 1 / 1 '). An important implication 
is the crossing of 8 curves at 4> c for different system sizes, 
provided the system size is large enough. Similar inter- 
sections are seen in the case of the Binder cumulanlj 1 ^ at 
the Curie temperature T c in the context of ferromagnetic 
transitions, and the crossing phenomena in both the for- 
mer and the latter represent a robust avenue for locating 
second order phase transitions. 

Fig. [T] shows 8 curves for different system sizes (i.e. 
L = 6.0, L = 9.0, and L = 13.5) with a well defined cross- 
ing at a common point. To determine 4> c in a systematic 
fashion, we calculate the order parameter for pairs of sys- 
tem sizes (e.g. {6, 9} and {8, 12}) with the larger member 
of each pair 1.5 times the size of the smaller. Although in 
the thermodynamic limit crossings would occur exactly 
at 4> c , the locations of intersections for finite size sys- 
tems [e.g. <j>c{L) differ slightly from <f) c due to finite size 
effects, and must be considered pscudocritical points 14 . 
The pseudo-critical excluded volume <p c (L) for each pair 
are recorded in table U and also appear in the graph in 
Fig. [21 with the system size reciprocal on the ab- 
scissa. For small system sizes, 4> C (L) rises, eventually 
saturating at 0.0317 ± 0.0004, the extrapolated result for 
the percolation threshold. 

A complementary finite size scaling technique for cal- 
culating the percolation threshold (j) c , based on the data 
collapse phenomenon in the critical region, simultane- 
ously takes into consideration results for a broad range 
of system sizes and also yields the correlation length ex- 
ponent v. 

Data collapse graphs may be prepared by plotting on 
the vertical axis the quantity 8 with L l / V \4> — <f> c \ (with 
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FIG. 2: (Color online) Pseudocritical excluded volumes (with 
calculated error bars) plotted with respect to the reciprocal 
of the mean system size L for each pair 



an optimal collapse for quantitatively correct values for 
v and (f> c ) on the abscissa for various systems sizes. The 
coincidence of data on the same curve is a manifestation 
of the existence of a universal singular scaling function 
for S as well as an indication that v and <p c are specified 
correctly. 

To use the data collapse as a quantitative tool, one 
notes that the "v" shaped curve in Fig. [3] is essentially 
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FIG. 3: (Color online) Data collapse with S data shown for 
various system sizes. The solid red line is a polynomial curve 
determined from a nonlinear least squares fit, and symbols 
represent data from Monte Carlo simulations. 
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FIG. 4: (Color online) \ plot shown with respect to the cor- 
relation length critical exponent v and the excluded volume 
<j)c with global minimum at v — 0.92 and <j> c = 0.0317. 



linear with a small degree of curvature. To accommodate 
this quasi-linear variation in a perturbative fashion, we 
use a form f(x) = Aq + A\x + • • • + A^x 4 " (the results are 
not impacted by the inclusion of the quartic term, and the 
series is truncated at fourth order) with x — L x l v ((j)—(f> c ). 
To assess the degree to which the points from Monte 
Carlo calculations lie on the same curve, we use 



X{<t>c,v, A , ■ ■ ■ , A 4 ) 



\ 



f(xi) - 5 t 



(1) 



The quantity \ provides an overall measure of the devia- 
tion of Monte Carlo data from the scaling curve. Hence, 
optimizing the data collapse entails minimizing \ with re- 
spect to cj> c , v, and the Ai parameters with a two tiered 
fit. First, for a particular choice of v and </> c , the param- 
eters Ai are calculated in a linear least squares fit. The 
appropriate values of the critical indices <j) c and v are 
then determined by a stochastically driven least squares 
fit, where a sequence of attempts are made to introduce 
small random perturbations to 4> c and is, with the shifts 
accepted only if the sum of squares term \ is thereby 
lowered. Fig. 0] shows x with respect to v and <J) C . The 
global minimum for <fi — 0.0317 and v — 0.92 is, within 
the bounds of Monte Carlo statistical error, consistent 
with the calculation of (j> c based on crossings of the 5 
quantity. Moreover, the calculated value of v is in accord 
with results gleaned from high precision Monte Carlo cal- 
culations for discrete 3D systems (e.g. 0.875 [1] for dis- 
crete percolation models on 3D lattices with cubic sym- 
metry^) . 



IV. CONCLUSIONS 

We have used a large-scale Monte Carlo calculation 
based on a dynamical exploration of void spaces to de- 
termine the critical properties (i.e. finding <ft c = 0.0317± 
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0.0004 and v = 0.92 ± 0.05) for the percolation through 
voids surrounding randomly placed impermeable spheres. 
Two complementary finite size scaling analyses yield the 
same percolation threshold <f> c . Our dynamically based 
calculation is distinct in that it makes no approximation 
for the structure of void spaces, while also incorporating 
finite size effects which are exploited to calculate criti- 
cal indices. The result for the correlation length critical 
exponent v is compatible with the critical behavior of 
discrete 3D counterparts. 

V. APPENDIX: SAMPLING THE POISSONIAN 
DISTRIBUTION 

To sample realizations of disorder from the appropriate 
statistical distribution, it is important to generate ran- 
dom configurations of spherical inclusions in an unbiased 
way; since we operate in the grand canonical ensemble, 
the number N of spheres in the simulation volume must 
in general vary from one sample to the next due to statis- 
tical fluctuations. The integer value closest to the mean 
occupancy iV av = pL 3 is a convenient initial choice, and 
random variations in the number of spheres in the system 
are taken into consideration with a sequence of stochas- 
tically driven attempts either to raise or lower N. The 
latter are part of an importance sampling scheme simi- 
lar to that used to derive the Metropolis criterion^! at 
the heart of Monte Carlo simulations that sample the 
Boltzmann distribution in the calculation of equilibrium 
thermodynamic variables. 

To determine the probability of having N spheres in 
the simulation volume v = L 3 , we divide v into M sub- 
volumes of equal size where Av = v/M. For large M the 
likelihood of multiple occupancy in any of the subvolumcs 



is very small relative to the chance of having one or zero 
sites in a subdivision; in the small Av limit, the single 
occupancy probability is pAv, with (1 — pAv) being the 
complementary likelihood of null occupancy. Hence, the 
probability that the entire system is devoid of hopping 
sites is (1 — pv/M) M , which becomes e~ pv for M —> oo. 
For single occupancy, adopting a prefactor M to take 
into consideration that the site may reside in any of the 
M subvolumes, yields M(pv/M)(l - pv/M) M - x , which 
becomes pve~ pv in the Aw — ► limit. Similar logic gives 
P{N) = e~ pv (pv) N /N\ for the general case of exactly N 
spheres in the simulation volume, where AH compensates 
for multiple occupancy. 

To generate a disorder realization, a succession of 
attempts (a number of moves in the vicinity of iV av 
is sufficient to achieve ergodicity) is made to raise or 
lower the occupancy number N, where the choice to in- 
crease or decrease N is randomly determined. For in- 
crements from N to N + 1, the change is accepted if 
X r < r+ = p{N + l)/p(N) = pv/(N + 1), where X r 
is a random number sampled uniformly from the inter- 
val [0,1]. Similarly, decreasing N to N — 1 occurs if 
X r < r_ = p(N - l)/p(N) = N/pv. Finally with N 
properly sampled, three Cartesian coordinates for each 
sphere are chosen independently (and at random with 
uniform probability density) from the interval [0,L]. 
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